####################################
#Shared Space and Civic Engagement #
# Replication Code #
## June 2024
####################################

rm(list = ls())

## set working directory to appropriate file path
#setwd("[your location path]/Replication materials")

## load packages
library(psych) 
library(ggplot2)
library(gridExtra)
library(merTools)
library(stargazer)
library(tidyverse)
library(DescTools)
library(sf)
library(spatstat)


#############################################
## Block parties on R's block
#############################################

## load in data set
data <- read.table("sharedspacedatav2.csv", head = TRUE, sep=",") ## 99,371 obs. of 23 vars

## 99,371 voters
unique(data$geoid_new) ## 381 census tracts

ls(data)

#############################################
#############################################
## Main paper
## ## Note: Figures 1 and 2 are conceptual (do not use data)
#############################################
#############################################


############################################
## Table 1: Association between party on 
## respondent’s block in 2012 and turnout (full sample)
#############################################

fit.basic2012 <- glmer(vote2012 ~ party2012 + 
                         (1 | geoid_new), family = binomial("logit"), 
                       control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.basic2012) 



fit.controls2012 <- glmer(vote2012 ~ party2012 + age + female + 
                            (1 | geoid_new), family = binomial("logit"), 
                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.controls2012) 


fit.full2012 <- glmer(vote2012 ~ party2012 + age + female +  
                        percrelchil + multi_unit_bin + med_rent_threehund +
                        perblackbin +  perforeignbornbin + dist_pol_mi + popdens_bin + perbachelorsbin + 
                        (1 | geoid_new), family = binomial("logit"), 
                      control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.full2012) 


## print to file: Table 1
stargazer(fit.basic2012, fit.controls2012, fit.full2012,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "html", 
          out = "Table1.html")

## Print in console: Table 1
stargazer(fit.basic2012, fit.controls2012, fit.full2012,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "text")

## Note: Census tract values drawn from model fit "groups": 381, 377, 372


############################################
## Figure 3:Predicted probability of turnout 
## by parties on a respondent’s block in 2012
#############################################

#############################################
## Prepare predicted probabilities for figure
#############################################

## at least one party in 2012
newdata <- with(data, data.frame(vote2012, party2012, geoid_new)) 


newdata[, 2] <- 1
head(newdata)

party <- predict(fit.basic2012, newdata, type="response",  allow.new.levels=TRUE)
head(party)
party.pred.naive <- mean(party,  na.rm = T)
party.pred.naive## 0.682


## no party in 2012
newdata <- with(data, data.frame(vote2012, party2012, geoid_new))

newdata[, 2] <- 0
head(newdata)

noparty <- predict(fit.basic2012, newdata, type="response",  allow.new.levels=TRUE)
noparty.pred.naive <- mean(noparty, na.rm = T)
noparty.pred.naive ## 0.669


## ind model

## at least one party in 2012
newdata <- with(data, data.frame(vote2012, party2012, geoid_new, age, female))

newdata[, 2] <- 1
newdata[, 4] <- median(data$age, na.rm = T)
newdata[, 5] <- median(data$female, na.rm = T)
head(newdata)

party <- predict(fit.controls2012, newdata, type="response", allow.new.levels=TRUE)
party.pred.ind <- mean(party, na.rm = T)
party.pred.ind## 0.645

## no party in 2012
newdata <- with(data, data.frame(vote2012, party2012, geoid_new, age, female))

newdata[, 2] <- 0
newdata[, 4] <- median(data$age, na.rm = T)
newdata[, 5] <- median(data$female, na.rm = T)
head(newdata)

noparty <- predict(fit.controls2012, newdata, type="response", allow.new.levels=TRUE)
noparty.pred.ind <- mean(noparty, na.rm = T)
noparty.pred.ind ## 0.626


## full model

## at least one party in 2012
newdata <- with(data, data.frame(vote2012, party2012, geoid_new, age, female, percrelchil, 
                                 multi_unit_bin, med_rent_threehund, ## tract predictors
                                 perblackbin,  perforeignbornbin, dist_pol_mi, popdens_bin, perbachelorsbin))

newdata[, 2] <- 1
newdata[, 4] <- median(data$age, na.rm = T)
newdata[, 5] <- median(data$female, na.rm = T)
newdata[, 6] <- median(data$percrelchil, na.rm = T)
newdata[, 7] <- median(data$multi_unit_bin, na.rm = T)
newdata[, 8] <- median(data$med_rent_threehund, na.rm = T)
newdata[, 9] <- median(data$perblackbin, na.rm = T)
newdata[, 10] <- median(data$perforeignbornbin, na.rm = T)
newdata[, 11] <- median(data$dist_pol_mi, na.rm = T)
newdata[, 12] <- median(data$popdens_bin, na.rm = T)
newdata[, 13] <- median(data$perbachelorsbin, na.rm = T)

head(newdata)

party <- predict(fit.full2012, newdata, type="response", allow.new.levels=TRUE)
party.pred.full <- mean(party, na.rm = T)
party.pred.full ## 0.689

## no party in 2012
newdata <- with(data, data.frame(vote2012, party2012, geoid_new, age, female, percrelchil, 
                                 multi_unit_bin, med_rent_threehund, ## tract predictors
                                 perblackbin,  perforeignbornbin, dist_pol_mi, popdens_bin, perbachelorsbin))


newdata[, 2] <- 0
newdata[, 4] <- median(data$age, na.rm = T)
newdata[, 5] <- median(data$female, na.rm = T)
newdata[, 6] <- median(data$percrelchil, na.rm = T)
newdata[, 7] <- median(data$multi_unit_bin, na.rm = T)
newdata[, 8] <- median(data$med_rent_threehund, na.rm = T)
newdata[, 9] <- median(data$perblackbin, na.rm = T)
newdata[, 10] <- median(data$perforeignbornbin, na.rm = T)
newdata[, 11] <- median(data$dist_pol_mi, na.rm = T)
newdata[, 12] <- median(data$popdens_bin, na.rm = T)
newdata[, 13] <- median(data$perbachelorsbin, na.rm = T)

head(newdata)

noparty <- predict(fit.full2012, newdata, type="response", allow.new.levels=TRUE)
noparty.pred.full <- mean(noparty, na.rm = T)

noparty.pred.full ## 0.674


## Comparison: Predicted probability of turnout by median rent
## high median rent
newdata <- with(data, data.frame(vote2012, party2012, geoid_new, age, female, percrelchil, 
                                 multi_unit_bin, med_rent_threehund, ## tract predictors
                                 perblackbin,  perforeignbornbin, dist_pol_mi, popdens_bin, perbachelorsbin))
newdata[, 2] <- median(data$party2012, na.rm = T)
newdata[, 4] <- median(data$age, na.rm = T)
newdata[, 5] <- median(data$female, na.rm = T)
newdata[, 6] <- median(data$percrelchil, na.rm = T)
newdata[, 7] <- median(data$multi_unit_bin, na.rm = T)
newdata[, 8] <- quantile(data$med_rent_threehund, 0.75, na.rm = T)
newdata[, 9] <- median(data$perblackbin, na.rm = T)
newdata[, 10] <- median(data$perforeignbornbin, na.rm = T)
newdata[, 11] <- median(data$dist_pol_mi, na.rm = T)
newdata[, 12] <- median(data$popdens_bin, na.rm = T)
newdata[, 13] <- median(data$perbachelorsbin, na.rm = T)


head(newdata)

high <- predict(fit.full2012, newdata, type="response", allow.new.levels=TRUE)
high.pred <- mean(high, na.rm = T)
high.pred## .695

## low median rent
newdata <- with(data, data.frame(vote2012, party2012, geoid_new, age, female, percrelchil, 
                                 multi_unit_bin, med_rent_threehund, ## tract predictors
                                 perblackbin,  perforeignbornbin, dist_pol_mi, popdens_bin, perbachelorsbin))
newdata[, 2] <- median(data$party2012, na.rm = T)
newdata[, 4] <- median(data$age, na.rm = T)
newdata[, 5] <- median(data$female, na.rm = T)
newdata[, 6] <- median(data$percrelchil,  na.rm = T)
newdata[, 7] <- median(data$multi_unit_bin, na.rm = T)
newdata[, 8] <- quantile(data$med_rent_threehund, 0.25, na.rm = T)
newdata[, 9] <- median(data$perblackbin, na.rm = T)
newdata[, 10] <- median(data$perforeignbornbin, na.rm = T)
newdata[, 11] <- median(data$dist_pol_mi, na.rm = T)
newdata[, 12] <- median(data$popdens_bin, na.rm = T)
newdata[, 13] <- median(data$perbachelorsbin, na.rm = T)


head(newdata)

low <- predict(fit.full2012, newdata, type="response", allow.new.levels=TRUE)
low.pred <- mean(low, na.rm = T)
low.pred## .685

############################################
## dataframe for plot
#############################################

data.test <- data.frame(
  Model = c("Naive model",
            "Naive model",
            "Individual model",
            "Individual model",
            "Full model",
            "Full model"),
  Party = c("At least one party", "No party"),
  effect = c(party.pred.naive, noparty.pred.naive,
             party.pred.ind, noparty.pred.ind,
             party.pred.full, noparty.pred.full)) 


data.test$Party <- factor(data.test$Party, levels =c("No party",
                                                     "At least one party"))
data.test


## Figure 3: Print in console
p1<- ggplot(data.test, aes(x=Party, y=effect, shape = Model, group = Model)) +
  geom_point(stat="identity", size = 3) +
  geom_line(linetype = "dashed") +
  ylab("Predicted probability of turnout") +
  guides(fill="none") +
  theme_bw() +
  theme(text = element_text(size=12) ,axis.title.x=element_blank()) +
  scale_y_continuous(limits=c(0.6, 0.7)) 
p1

## Figure 3: Print PDF
pdf("figure3.pdf", width = 6, height = 3)
grid.arrange(p1)
dev.off()

#############################################
## Figure 4: Distribution of voters without a block party within 152 meters 
## in 2011 by presence/absence of a block party within 152 meters in 2012 
## (quasi-experimental test one)
#############################################
dat_vot<- st_read("full100_party_wgeo_andlatlon_reduced.geojson") #full data with yes/no indicators for various party permutations and spatial location

parties<-st_read("blkparty_points1_copy.geojson") #spatial location of block parties

phl_lim<-st_read( "http://data.phl.opendata.arcgis.com/datasets/405ec3da942d4e20869d4e1449a2be48_0.geojson") #philadelphia spatial boundary

#st_transform(phl_lim) = 2272 #setting projection for Philadelphia boundary

#setting the voter location, party location, and Philadelphia boundary to same projection
dat_vot_prj<-st_transform(dat_vot, 32618)
st_crs(dat_vot_prj)
parties_prj<-st_transform(parties, 32618)

phl_lim_prj<-st_transform(phl_lim, 32618)

plot(dat_vot_prj[1]) #visual check to make sure voters' locations are correct -- note that some voters appear to be outside Philly city limits

dat_vot2<-st_intersection(dat_vot_prj, phl_lim_prj) #subsetting voter data to residential locations only within boundaries of Philadelphia



p2 <- ggplot(dat_vot2[dat_vot2$party2011==0,]) + geom_sf(aes(color = as.factor(party12n11))) +  geom_sf(data = phl_lim_prj, fill=NA, color = "black") +
  scale_color_grey(name = "Voter by Treatment", breaks = c("0", "1"), 
                       labels = c("Control: No Party in 2011, No Party in 2012", "Treatement: No Party in 2011, Party in 2012"),   
                   start = 0.8, end = 0.2) +
  theme(axis.text = element_blank(), axis.ticks = element_blank()) 

p2


## Figure 4: Print PDF
pdf("figure4.pdf", width = 12, height = 8)
grid.arrange(p2)
dev.off()
#############################################
## Figure 5: Distribution of voters with a block party within 152 meters 
## in 2011 by presence/absence of a block party within 152 meters in 2012 
## (quasi-experimental test two)
#############################################

#note this figure builds on data processing for figure 4

p3 <- ggplot(dat_vot2[dat_vot2$party2011==1,]) + geom_sf(aes(color = as.factor(party11n12))) +  geom_sf(data = phl_lim_prj, fill=NA, color = "black") +
  scale_color_grey(name = "Voter by Treatment", breaks = c("0", "1"), 
                       labels = c("Control: Party in 2011, Party in 2012", "Treatment: Party in 2011, No Party in 2012"),
                   start = 0.8, end = 0.2) +
  theme(axis.text = element_blank(), axis.ticks = element_blank())

p3

## Figure 5: Print PDF
pdf("figure5.pdf", width = 12, height = 8)
grid.arrange(p3)
dev.off()

#############################################
## Table 2: Comparing demographic characteristics 
## across treated and control units
## Note: fill values into table
#############################################

## Subset 1: No party on R’s block in 2011
data.noparty2011 <- data[data$party2011 == 0, ] ## 26,951 obs.
table(data.noparty2011$vote2012)
prop.table(table(data.noparty2011$vote2012))

## compare treated and control units

## Quasi-experimental test 1

## perc. related child
t.test(data.noparty2011$percrelchil[data.noparty2011$party2012 == 1], 
       data.noparty2011$percrelchil[data.noparty2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) 

## perc. multi-unit
t.test(data.noparty2011$multi_unit[data.noparty2011$party2012 == 1], 
       data.noparty2011$multi_unit[data.noparty2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) 

## median rent
t.test(data.noparty2011$med_rent_gE[data.noparty2011$party2012 == 1], 
       data.noparty2011$med_rent_gE[data.noparty2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) 

## perc. African American
t.test(data.noparty2011$perblack[data.noparty2011$party2012 == 1], 
       data.noparty2011$perblack[data.noparty2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) 

## perc. foreign born
t.test(data.noparty2011$perforeignborn[data.noparty2011$party2012 == 1], 
       data.noparty2011$perforeignborn[data.noparty2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) 

## perc. bachelor's degree
t.test(data.noparty2011$perbachelors[data.noparty2011$party2012 == 1], 
       data.noparty2011$perbachelors[data.noparty2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) 

## distance from polling place
t.test(data.noparty2011$dist_pol_mi[data.noparty2011$party2012 == 1], 
       data.noparty2011$dist_pol_mi[data.noparty2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) 

## population density
t.test(data.noparty2011$popden[data.noparty2011$party2012 == 1], 
       data.noparty2011$popden[data.noparty2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) 

## Subset 2: At least one party on R’s block in 2011

data.party2011 <- data[data$party2011 == 1, ] ## 72,420 obs
table(data.party2011$vote2012)
prop.table(table(data.party2011$vote2012))

## create indicator for no party in 2012
data.party2011$noparty2012 <- 0
data.party2011$noparty2012[data.party2011$party2012 ==0] <- 1
table(data.party2011$noparty2012)

## compare treated and control units

## Quasi-experimental test 2

## perc. related child
t.test(data.party2011$percrelchil[data.party2011$party2012 == 1], 
       data.party2011$percrelchil[data.party2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) #sig

## perc. multi-unit
t.test(data.party2011$multi_unit[data.party2011$party2012 == 1], 
       data.party2011$multi_unit[data.party2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) #sig

## median rent
t.test(data.party2011$med_rent_gE[data.party2011$party2012 == 1], 
       data.party2011$med_rent_gE[data.party2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) #sig

## perc. African American
t.test(data.party2011$perblack[data.party2011$party2012 == 1], 
       data.party2011$perblack[data.party2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) #sig

## perc. foreign born
t.test(data.party2011$perforeignborn[data.party2011$party2012 == 1], 
       data.party2011$perforeignborn[data.party2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) #sig

## perc. bachelor's
t.test(data.party2011$perbachelors[data.party2011$party2012 == 1], 
       data.party2011$perbachelors[data.party2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) #sig

## distance from polling place
t.test(data.party2011$dist_pol_mi[data.party2011$party2012 == 1], 
       data.party2011$dist_pol_mi[data.party2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) #sig

## population density
t.test(data.party2011$popden[data.party2011$party2012 == 1], 
       data.party2011$popden[data.party2011$party2012 == 0],
       alternative = c("two.sided"),  var.equal = TRUE, conf.level = 0.95) #sig

#############################################
## Figure 6: Effects of block parties on turnout in 2012 (quasi-experiments)
#############################################

############################################
## Quasi-experiment 1
#############################################

fit.basic2012.ne1 <- glmer(vote2012 ~ party2012 + 
                         (1 | geoid_new), family = binomial("logit"), 
                       control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.noparty2011)

summary(fit.basic2012.ne1) 



fit.controls2012.ne1 <- glmer(vote2012 ~ party2012 + age + female + 
                            (1 | geoid_new), family = binomial("logit"), 
                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.noparty2011)

summary(fit.controls2012.ne1) 


fit.full2012.ne1 <- glmer(vote2012 ~ party2012 + age + female + 
                        percrelchil + multi_unit_bin + med_rent_threehund +
                        perblackbin +  perforeignbornbin + dist_pol_mi +  popdens_bin + perbachelorsbin + 
                        (1 | geoid_new), family = binomial("logit"), 
                      control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.noparty2011)

summary(fit.full2012.ne1) 

#############################################
## QE1: Predicted probabilities for plot
#############################################

## at least one party in 2012
newdata <- with(data.noparty2011, data.frame(vote2012, party2012, geoid_new))


newdata[, 2] <- 1
head(newdata)

party <- predict(fit.basic2012.ne1, newdata, type="response",  allow.new.levels=TRUE)
head(party)
party.pred.naive <- mean(party,  na.rm = T)
party.pred.naive ## 0.695


## no party in 2012
newdata <- with(data.noparty2011, data.frame(vote2012, party2012, geoid_new))

newdata[, 2] <- 0
head(newdata)

noparty <- predict(fit.basic2012.ne1, newdata, type="response",  allow.new.levels=TRUE)
noparty.pred.naive <- mean(noparty, na.rm = T)
noparty.pred.naive ## 0.672


## ind model

## at least one party in 2012
newdata <- with(data.noparty2011, data.frame(vote2012, party2012, geoid_new, age, female))

newdata[, 2] <- 1
newdata[, 4] <- median(data.noparty2011$age, na.rm = T)
newdata[, 5] <- median(data.noparty2011$female, na.rm = T)
head(newdata)

party <- predict(fit.controls2012.ne1, newdata, type="response", allow.new.levels=TRUE)
party.pred.ind <- mean(party, na.rm = T)
party.pred.ind ## 0.670

## no party in 2012
newdata <- with(data.noparty2011, data.frame(vote2012, party2012, geoid_new, age, female))

newdata[, 2] <- 0
newdata[, 4] <- median(data.noparty2011$age, na.rm = T)
newdata[, 5] <- median(data.noparty2011$female, na.rm = T)
head(newdata)

noparty <- predict(fit.controls2012.ne1, newdata, type="response", allow.new.levels=TRUE)
noparty.pred.ind <- mean(noparty, na.rm = T)
noparty.pred.ind ## 0.644


## full model

## at least one party in 2012
newdata <- with(data.noparty2011, data.frame(vote2012, party2012, geoid_new, age, female, percrelchil, 
                                             multi_unit_bin, med_rent_threehund, ## tract predictors
                                             perblackbin,  perforeignbornbin, dist_pol_mi, popdens_bin, perbachelorsbin))

newdata[, 2] <- 1
newdata[, 4] <- median(data.noparty2011$age, na.rm = T)
newdata[, 5] <- median(data.noparty2011$female, na.rm = T)
newdata[, 6] <- median(data.noparty2011$percrelchil, na.rm = T)
newdata[, 7] <- median(data.noparty2011$multi_unit_bin, na.rm = T)
newdata[, 8] <- median(data.noparty2011$med_rent_threehund, na.rm = T)
newdata[, 9] <- median(data.noparty2011$perblackbin, na.rm = T)
newdata[, 10] <- median(data.noparty2011$perforeignbornbin, na.rm = T)
newdata[, 11] <- median(data.noparty2011$dist_pol_mi, na.rm = T)
newdata[, 12] <- median(data.noparty2011$popdens_bin, na.rm = T)
newdata[, 13] <- median(data.noparty2011$perbachelorsbin, na.rm = T)

head(newdata)

party <- predict(fit.full2012.ne1, newdata, type="response", allow.new.levels=TRUE)
party.pred.full <- mean(party, na.rm = T)
party.pred.full ## 0.641

## no party in 2012
newdata <- with(data.noparty2011, data.frame(vote2012, party2012, geoid_new, age, female, percrelchil, 
                                             multi_unit_bin, med_rent_threehund, ## tract predictors
                                             perblackbin,  perforeignbornbin, dist_pol_mi, popdens_bin, perbachelorsbin))


newdata[, 2] <- 0
newdata[, 4] <- median(data.noparty2011$age, na.rm = T)
newdata[, 5] <- median(data.noparty2011$female, na.rm = T)
newdata[, 6] <- median(data.noparty2011$percrelchil, na.rm = T)
newdata[, 7] <- median(data.noparty2011$multi_unit_bin, na.rm = T)
newdata[, 8] <- median(data.noparty2011$med_rent_threehund, na.rm = T)
newdata[, 9] <- median(data.noparty2011$perblackbin, na.rm = T)
newdata[, 10] <- median(data.noparty2011$perforeignbornbin, na.rm = T)
newdata[, 11] <- median(data.noparty2011$dist_pol_mi, na.rm = T)
newdata[, 12] <- median(data.noparty2011$popdens_bin, na.rm = T)
newdata[, 13] <- median(data.noparty2011$perbachelorsbin, na.rm = T)

head(newdata)

noparty <- predict(fit.full2012.ne1, newdata, type="response", allow.new.levels=TRUE)
noparty.pred.full <- mean(noparty, na.rm = T)
noparty.pred.full ## 0.613

############################################
## plot: Panel 1 (Quasi-Experiment 1)
#############################################

data.test <- data.frame(
  Model = c("Naive model",
            "Naive model",
            "Individual model",
            "Individual model",
            "Full model",
            "Full model"),
  Party = c("At least one party", "No party"),
  effect = c(party.pred.naive, noparty.pred.naive,
             party.pred.ind, noparty.pred.ind,
             party.pred.full, noparty.pred.full)) 


data.test$Party <- factor(data.test$Party, levels =c("No party",
                                                     "At least one party"))
data.test


p4<- ggplot(data.test, aes(x=Party, y=effect, shape = Model, group = Model)) +
  geom_point(stat="identity", size = 3) +
  geom_line(linetype = "dashed") +
  ylab("Predicted probability of turnout") +
  guides(fill="none") +
  theme_bw() +
  theme(text = element_text(size=12) ,axis.title.x=element_blank()) +
  scale_y_continuous(limits=c(0.6, 0.7)) +
  ggtitle("Quasi-experimental test one")
p4


############################################
## Quasi-experiment 2
#############################################
## Note: created this indicator for Table 2
## create indicator for no party in 2012
## data.party2011$noparty2012 <- 0
## data.party2011$noparty2012[data.party2011$party2012 ==0] <- 1
## table(data.party2011$noparty2012)

fit.basic2012.ne2 <- glmer(vote2012 ~ noparty2012 + 
                         (1 | geoid_new), family = binomial("logit"), 
                       control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.party2011)

summary(fit.basic2012.ne2) 



fit.controls2012.ne2 <- glmer(vote2012 ~ noparty2012 + age + female +  
                            (1 | geoid_new), family = binomial("logit"), 
                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.party2011)

summary(fit.controls2012.ne2) 


fit.full2012.ne2 <- glmer(vote2012 ~ noparty2012 + age + female +  
                        percrelchil + multi_unit_bin + med_rent_threehund +
                        perblackbin +  perforeignbornbin + dist_pol_mi +  popdens_bin + perbachelorsbin +
                        (1 | geoid_new), family = binomial("logit"), 
                      control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.party2011)

summary(fit.full2012.ne2) 

#############################################
## QE2: Predicted probabilities for plot
#############################################


## no party in 2012
newdata <- with(data.party2011, data.frame(vote2012, noparty2012, geoid_new)) 


newdata[, 2] <- 1
head(newdata)

party <- predict(fit.basic2012.ne2, newdata, type="response",  allow.new.levels=TRUE)
head(party)
party.pred.naive <- mean(party,  na.rm = T)
party.pred.naive## 0.670


## at least one party in 2012
newdata <- with(data.party2011, data.frame(vote2012, noparty2012, geoid_new))

newdata[, 2] <- 0
head(newdata)

noparty <- predict(fit.basic2012.ne2, newdata, type="response",  allow.new.levels=TRUE)
noparty.pred.naive  <- mean(noparty, na.rm = T)
noparty.pred.naive ## 0.680


## ind model

## no party in 2012
newdata <- with(data.party2011, data.frame(vote2012, noparty2012, geoid_new, age, female))

newdata[, 2] <- 1
newdata[, 4] <- median(data.party2011$age, na.rm = T)
newdata[, 5] <- median(data.party2011$female, na.rm = T)
head(newdata)

party <- predict(fit.controls2012.ne2, newdata, type="response", allow.new.levels=TRUE)
party.pred.ind  <- mean(party, na.rm = T)
party.pred.ind # 0.626

## at least one party in 2012
newdata <- with(data.party2011, data.frame(vote2012, noparty2012, geoid_new, age, female))

newdata[, 2] <- 0
newdata[, 4] <- median(data.party2011$age, na.rm = T)
newdata[, 5] <- median(data.party2011$female, na.rm = T)
head(newdata)

noparty <- predict(fit.controls2012.ne2, newdata, type="response", allow.new.levels=TRUE)
noparty.pred.ind  <- mean(noparty, na.rm = T)
noparty.pred.ind ## 0.638


## full model

## no party in 2012
newdata <- with(data.party2011, data.frame(vote2012, noparty2012, geoid_new, age, female, percrelchil, 
                                           multi_unit_bin, med_rent_threehund, ## tract predictors
                                           perblackbin,  perforeignbornbin, dist_pol_mi, popdens_bin, perbachelorsbin ))

newdata[, 2] <- 1
newdata[, 4] <- median(data.party2011$age, na.rm = T)
newdata[, 5] <- median(data.party2011$female, na.rm = T)
newdata[, 6] <- median(data.party2011$percrelchil, na.rm = T)
newdata[, 7] <- median(data.party2011$multi_unit_bin, na.rm = T)
newdata[, 8] <- median(data.party2011$med_rent_threehund, na.rm = T)
newdata[, 9] <- median(data.party2011$perblackbin, na.rm = T)
newdata[, 10] <- median(data.party2011$perforeignbornbin, na.rm = T)
newdata[, 11] <- median(data.party2011$dist_pol_mi, na.rm = T)
newdata[, 12] <- median(data.party2011$popdens_bin, na.rm = T)
newdata[, 13] <- median(data.party2011$perbachelorsbin, na.rm = T)

head(newdata)

party <- predict(fit.full2012.ne2, newdata, type="response", allow.new.levels=TRUE)
party.pred.full <- mean(party, na.rm = T)
party.pred.full ## 0.675

## no party in 2012
newdata <- with(data.party2011, data.frame(vote2012, noparty2012, geoid_new, age, female, percrelchil, 
                                           multi_unit_bin, med_rent_threehund, ## tract predictors
                                           perblackbin,  perforeignbornbin, dist_pol_mi, popdens_bin, perbachelorsbin))


newdata[, 2] <- 0
newdata[, 4] <- median(data.party2011$age, na.rm = T)
newdata[, 5] <- median(data.party2011$female, na.rm = T)
newdata[, 6] <- median(data.party2011$percrelchil, na.rm = T)
newdata[, 7] <- median(data.party2011$multi_unit_bin, na.rm = T)
newdata[, 8] <- median(data.party2011$med_rent_threehund, na.rm = T)
newdata[, 9] <- median(data.party2011$perblackbin, na.rm = T)
newdata[, 10] <- median(data.party2011$perforeignbornbin, na.rm = T)
newdata[, 11] <- median(data.party2011$dist_pol_mi, na.rm = T)
newdata[, 12] <- median(data.party2011$popdens_bin, na.rm = T)
newdata[, 13] <- median(data.party2011$perbachelorsbin, na.rm = T)

head(newdata)

noparty <- predict(fit.full2012.ne2, newdata, type="response", allow.new.levels=TRUE)
noparty.pred.full <- mean(noparty, na.rm = T)
noparty.pred.full ## 0.685

############################################
## plot: Panel 2 (Quasi-Experiment 2)
#############################################

data.test <- data.frame(
  Model = c("Naive model",
            "Naive model",
            "Individual model",
            "Individual model",
            "Full model",
            "Full model"),
  Party = c("No party", "At least one party"),
  effect = c(party.pred.naive, noparty.pred.naive,
             party.pred.ind, noparty.pred.ind,
             party.pred.full, noparty.pred.full)) 


data.test


p5<- ggplot(data.test, aes(x=Party, y=effect, shape = Model, group = Model)) +
  geom_point(stat="identity", size = 3) +
  geom_line(linetype = "dashed") +
  ylab("Predicted probability of turnout") +
  guides(fill="none") +
  theme_bw() +
  theme(text = element_text(size=12) ,axis.title.x=element_blank()) +
  scale_y_continuous(limits=c(0.6, 0.7)) +
  ggtitle("Quasi-experimental test two")
p5

grid.arrange(p4, p5, nrow =1)

## Figure 6: Print PDF
pdf("figure6.pdf", width = 12, height = 6)
grid.arrange(p4, p5, nrow =1)
dev.off()

#############################################
## Table 3: Association between party on respondent’s block 
## in 2012 and turnout by percentage of African American residents in census tract
#############################################

## indicator for high per black: cut-off at median value for subset
data$perblackmedian <- ifelse(data$perblack >= 32, 1, 0)
table(data$perblackmedian)

############################################
## Regression model
## percent black above median level
#############################################

fit.full2012.highperblack <- glmer(vote2012 ~ party2012 + age + female +  
                                     percrelchil + multi_unit_bin + med_rent_threehund + 
                                     perforeignbornbin + dist_pol_mi + popdens_bin + perbachelorsbin +  
                                     (1 | geoid_new), family = binomial("logit"), 
                                   control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data[data$perblackmedian == 1,])

summary(fit.full2012.highperblack) 

############################################
## Regression model
## percent black below median level
#############################################

fit.full2012.lowperblack <- glmer(vote2012 ~ party2012 + age + female +  
                                    percrelchil + multi_unit_bin + med_rent_threehund + 
                                    perforeignbornbin + dist_pol_mi + popdens_bin + perbachelorsbin +
                                    (1 | geoid_new), family = binomial("logit"), 
                                  control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data[data$perblackmedian == 0,])

summary(fit.full2012.lowperblack) 

############################################
## Regression model## interactive model
#############################################

fit.full2012.int <- glmer(vote2012 ~ party2012 + age + female + 
                            percrelchil + multi_unit_bin + med_rent_threehund +
                            perblackbin +  perforeignbornbin + dist_pol_mi + popdens_bin + perbachelorsbin +
                            perblackbin*party2012+
                            (1 | geoid_new), family = binomial("logit"), 
                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.full2012.int)

## Print to file: Table 3
stargazer(fit.full2012.highperblack, fit.full2012.lowperblack, fit.full2012.int,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelor's degree (binary)",
                               "At lest one party on R's block in 2012 x Percent black",
                               "Constant"
                               
          ), type = "html", 
          out = "Table3.html")

## Print to console: Table 3
stargazer(fit.full2012.highperblack, fit.full2012.lowperblack, fit.full2012.int,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelor's degree (binary)",
                               "At lest one party on R's block in 2012 x Percent black",
                               "Constant"
                               
          ), type = "text")

## Note: Census tract values drawn from model fit "groups": 184, 188, 372

#############################################
#############################################
## Appendix
## Note: Table A2 presents coding scheme (does not use data)
#############################################
#############################################

#############################################
## Table A1: Block party event type (2006-2016)
#############################################
dat <- read.csv("blkparty_points1_copy.csv", header = T) ## 105,610 obs. of block parties
ls(dat)
table(dat$event_type)

dat<- dat %>% mutate(even_recode = 
                       case_match(
                         event_type, 
                         "4TH OF JULY" ~ "Community",
                         "BIRTHDAY PARTY" ~ "Personal",
                         "COMMUNION" ~ "Community",
                         "EASTER EGG HUNT" ~ "Community",
                         "GRADUATION PARTY" ~ "Personal",
                         "LABOR DAY" ~ "Community",
                         "MOTHER DAY" ~ "Community",
                         "NEW YEARS EVE" ~ "Community",
                         "REPASS" ~ "Other",
                         "ARTS & CRAFTS SHOW" ~ "Community",
                         "STOP THE VIOLENCE CRUSADE" ~ "Community",
                         "CHRISTMAS PARTY" ~ "Community",
                         "COMMUNITY FUN DAY" ~ "Community",
                         "FATHER DAY" ~ "Community",
                         "HALLOWEEN PARTY" ~ "Community",
                         "MAY DAY" ~ "Community",
                         "NATIONAL NIGHT OUT" ~ "Community",
                         "PROM" ~ "Personal",
                         "SERENADES" ~ "Personal",
                         "WEDDING" ~ "Personal",
                         "BABY SHOWER" ~ "Personal",
                         "CHURCH SERVICE" ~ "Community",
                         "DEDICATION" ~ "Community",
                         "GENERAL" ~ "General",
                         "HEALTH FAIR" ~ "Community",
                         "MEMORIAL DAY" ~ "Community",
                         "NEW YEARS DAY" ~ "Community",
                         "RELIGIOUS EVENT" ~ "Community",
                         "SPRING FESTIVAL" ~ "Community",
                         "WEDDING RECEPTION" ~ "Personal"
                         
                       ))
Freq(dat$even_recode)

#############################################
## Table A3: Summary Statistics
#############################################

describe(data$multi_unit) ## perc. multi unit
describe(data$percrelchil) ## perc. related child
describe(data$med_rent_threehund) ## median rent
describe(data$popden) ## population density
describe(data$perblack) ## perc. African American
describe(data$perforeignborn) ## perc. foreign born
describe(data$perbachelors) ## perc. bachelors
describe(data$dist_pol_mi) ## distance from polling place
describe(data$age) ## age
describe(data$female) ## gender
describe(data$party2011) ## block party (2011)
describe(data$party2012) ## block party (2012)
describe(data$vote2012) ## voted in 2012

############################################
# Figure A1: Distributions of key variables
############################################

colnames(data)
dat_dist <- data |> select(multi_unit, percrelchil, med_rent_threehund, popden, perblack, perforeignborn, 
                           perbachelors, dist_pol_mi, age, female, party2011, party2011, party2012, vote2012) |>
  rename(
    "Percent multi-unit" = multi_unit, 
    "Percent related child" = percrelchil,
    "Median rent (unit = $300)" = med_rent_threehund,
    "Population density" = popden,
    "Percent African American" = perblack,
    "Percent foreign born" = perforeignborn,
    "Percent bachelors" = perbachelors,
    "Distance from nearest polling place" = dist_pol_mi,
    "Age" = age,
    "Female" = female, 
    "Block party (2011)" = party2011,
    "Block party (2012)" = party2012,
    "Vote" = vote2012
    )

dat_dist_plot <- dat_dist |> pivot_longer(cols = everything())

p6 <- ggplot(dat_dist_plot, aes(x = value)) + geom_histogram() + facet_wrap(~name, scales = "free")
p6

## Figure A1: Print PDF
pdf("figureA1.pdf", width = 10, height = 8)
grid.arrange(p6, nrow =1)
dev.off()

############################################
## Figure A2: Distribution of block parties in 2011 and 2012
############################################

#note this builds on data processing in Figure 4

parties_prj$year<- format(parties_prj$start_date,"%Y")
#head(parties$year)



#part_sub<- parties[parties$year==c("2011", "2012"),]
part_sub<-parties_prj[parties_prj$year=="2011" | parties_prj$year=="2012",]

part_sub$party<-1


p7 <- ggplot(data = part_sub) + geom_sf(aes(color = as.factor(party))) + geom_sf(data = phl_lim_prj, fill=NA, color = "black") + facet_wrap(.~ year) + scale_color_discrete(name = "", labels = " = 1 Party") + theme(axis.text = element_blank(), axis.ticks = element_blank())
p7

## Figure A2: Print PDF
pdf("figureA2.pdf", width = 12, height = 6)
grid.arrange(p7, nrow =1)
dev.off()

############################################
## Figure A3: Kernel density (Treatment: Party in 2011, No Party in 2012), Sigma = 500 Meters
############################################
#note: builds on data processing from Figure A2 and Figure 4

party11_ppp<- as.ppp(part_sub[part_sub$year=="2011",])
marks(party11_ppp)<-NULL
party12_ppp<- as.ppp(part_sub[part_sub$year=="2012",])

treat1<- as.ppp(dat_vot2[dat_vot2$party2011==1 & dat_vot2$party11n12==1,])
marks(treat1)<-NULL

treat2<- as.ppp(dat_vot2[dat_vot2$party2011==0 & dat_vot2$party12n11==1,])
marks(treat2)<-NULL

phl_win<- as.owin(phl_lim_prj)


Window(treat1)<- phl_win
plot(treat1)
treat1_kern<-density(treat1, sigma = 500)
plot(treat1_kern, main = NULL)


## Figure A3: Print PDF
pdf("figureA3.pdf", width = 10, height = 8)
plot(treat1_kern, main = NULL)
dev.off()

############################################
## Figure A4: Kernel density (Treatment: No party in 2011, Party in 2012), Sigma = 500 Meters
############################################
#note: builds on data processing from Figure A3

Window(treat2)<- phl_win
plot(treat2)
treat2_kern<-density(treat2, sigma = 500)
plot(treat2_kern, main = NULL)

## Figure A4: Print PDF
pdf("figureA4.pdf", width = 10, height = 8)
plot(treat2_kern, main = NULL)
dev.off()

############################################
## Figure A5: Simulation of complete spatial randomness (Treatment: Party in 2011, no Party in 2012)
############################################
#note: builds on data processing from Figure A3

l_treat1<- envelope(treat1, Lest, correction="border") #envelope will vary slightly based on random nature of the simulations
plot(l_treat1, main = NULL) 

## Figure A5: Print PDF
pdf("figureA5.pdf", width = 10, height = 8)
plot(l_treat1, main = NULL) 
dev.off()

############################################
## Figure A6: Simulation of complete spatial randomness (Treatment: Party in 2012, no Party in 2011)
############################################
#note: builds on data processing from Figure A4

l_treat2<- envelope(treat2, Lest, correction="border") #envelope will vary slightly based on random nature of the simulations
plot(l_treat2, main = NULL)

## Figure A6: Print PDF
pdf("figureA6.pdf", width = 10, height = 8)
plot(l_treat2, main = NULL) 
dev.off()

############################################
## Table A4: Distribution of turnout in 2012 across quasi-experimental tests
## Note: fill values into table
############################################
## Subset 1: No party on R’s block in 2011
data.noparty2011 <- data[data$party2011 == 0, ] ## 26951 obs.

## voting overall
table(data.noparty2011$vote2012)
prop.table(table(data.noparty2011$vote2012))

## treated and control Ns
table(data.noparty2011$party2012) ## 5696 of 26718

## voted among treated
table(data.noparty2011$vote2012[data.noparty2011$party2012 == 1])
prop.table(table(data.noparty2011$vote2012[data.noparty2011$party2012 == 1]))

## voting among control
table(data.noparty2011$vote2012[data.noparty2011$party2012 == 0])
prop.table(table(data.noparty2011$vote2012[data.noparty2011$party2012 == 0]))

## Subset 2: At least one party on R’s block in 2011
data.party2011 <- data[data$party2011 == 1, ] ## 72,420 obs.

## voting overall
table(data.party2011$vote2012)
prop.table(table(data.party2011$vote2012))

## treated and control Ns (note: in this case the treatment = 0)
table(data.party2011$party2012)

## create indicator for "no party in 2012" (treatment for subset 2 now = 1)
data.party2011$noparty2012 <- 0
data.party2011$noparty2012[data.party2011$party2012 ==0] <- 1
table(data.party2011$noparty2012)

## voted among treated
table(data.party2011$vote2012[data.party2011$noparty2012 == 1])
prop.table(table(data.party2011$vote2012[data.party2011$noparty2012 == 1]))

## voted among control
table(data.party2011$vote2012[data.party2011$noparty2012 == 0])
prop.table(table(data.party2011$vote2012[data.party2011$noparty2012 == 0]))

#############################################
## Table A5: Association between party on respondent’s block 
## in 2012 and turnout (Quasi-experimental test one)
#############################################

############################################
## Quasi-experiment 1
## Note: regression models from Figure 6
#############################################
## Note: code builds on commented models (from Figure 6)
#fit.basic2012.ne1 <- glmer(vote2012 ~ party2012 + 
#                             (1 | geoid_new), family = binomial("logit"), 
#                           control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.noparty2011)

summary(fit.basic2012.ne1) 



#fit.controls2012.ne1 <- glmer(vote2012 ~ party2012 + age + female + 
#                                (1 | geoid_new), family = binomial("logit"), 
#                              control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.noparty2011)

summary(fit.controls2012.ne1) 


#fit.full2012.ne1 <- glmer(vote2012 ~ party2012 + age + female + 
#                            percrelchil + multi_unit_bin + med_rent_threehund +
#                            perblackbin +  perforeignbornbin + dist_pol_mi +  popdens_bin + perbachelorsbin + 
#                            (1 | geoid_new), family = binomial("logit"), 
#                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.noparty2011)

summary(fit.full2012.ne1) 

## Print to file: Table A5
stargazer(fit.basic2012.ne1, fit.controls2012.ne1, fit.full2012.ne1,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "html", 
          out = "TableA5.html")

## Print to console: Table A5
stargazer(fit.basic2012.ne1, fit.controls2012.ne1, fit.full2012.ne1,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "text")

## Note: Census tract values drawn from model fit "groups": 317, 311, 307

#############################################
## Table A6: Association between no party on respondent’s block in 2012 
## and turnout (Quasi-experimental test two)
#############################################

############################################
## Quasi-experiment 2
## Note: regression models from Figure 6
#############################################

#fit.basic2012.ne2 <- glmer(vote2012 ~ noparty2012 + 
#                             (1 | geoid_new), family = binomial("logit"), 
#                           control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.party2011)

summary(fit.basic2012.ne2) 



#fit.controls2012.ne2 <- glmer(vote2012 ~ noparty2012 + age + female +  
#                                (1 | geoid_new), family = binomial("logit"), 
#                              control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.party2011)

summary(fit.controls2012.ne2) 


#fit.full2012.ne2 <- glmer(vote2012 ~ noparty2012 + age + female +  
#                            percrelchil + multi_unit_bin + med_rent_threehund +
#                            perblackbin +  perforeignbornbin + dist_pol_mi +  popdens_bin + perbachelorsbin +
#                            (1 | geoid_new), family = binomial("logit"), 
#                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.party2011)

summary(fit.full2012.ne2) 


## Print to file: Table A6
stargazer(fit.basic2012.ne2, fit.controls2012.ne2, fit.full2012.ne2,
          covariate.labels = c("No party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "html", 
          out = "TableA6.html")

## Print to console: Table A6
stargazer(fit.basic2012.ne2, fit.controls2012.ne2, fit.full2012.ne2,
          covariate.labels = c("No party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "text")

## Note: Census tract values drawn from model fit "groups": 349, 348, 347

############################################
## Table A7: Odds ratios for association between 
## party on respondent’s block in 2012 and turnout (full sample)
## Note: regression models from Table 1 and Figure 3
#############################################

#fit.basic2012 <- glmer(vote2012 ~ party2012 + 
#                         (1 | geoid_new), family = binomial("logit"), 
#                       control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.basic2012) 



#fit.controls2012 <- glmer(vote2012 ~ party2012 + age + female + 
#                            (1 | geoid_new), family = binomial("logit"), 
#                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.controls2012) 


#fit.full2012 <- glmer(vote2012 ~ party2012 + age + female +  
#                        percrelchil + multi_unit_bin + med_rent_threehund +
#                        perblackbin +  perforeignbornbin + dist_pol_mi + popdens_bin + perbachelorsbin + 
#                        (1 | geoid_new), family = binomial("logit"), 
#                      control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.full2012) 

## Note: removed standard errors, stars, and "Constant" from Table A7 in the appendix

## Print to file: Table A7
stargazer(fit.basic2012, fit.controls2012, fit.full2012,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "html", apply.coef = exp,
          out = "TableA7.html")

## Print in console: Table A7
stargazer(fit.basic2012, fit.controls2012, fit.full2012,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "text", apply.coef = exp)

############################################
## Table A8: Odds ratios for association between party on 
## respondent’s block in 2012 and turnout (Quasi-experimental test one)
############################################

############################################
## Quasi-experiment 1
## Note: regression models from Figure 6
#############################################

#fit.basic2012.ne1 <- glmer(vote2012 ~ party2012 + 
#                             (1 | geoid_new), family = binomial("logit"), 
#                           control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.noparty2011)

summary(fit.basic2012.ne1) 



#fit.controls2012.ne1 <- glmer(vote2012 ~ party2012 + age + female + 
#                                (1 | geoid_new), family = binomial("logit"), 
#                              control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.noparty2011)

summary(fit.controls2012.ne1) 


#fit.full2012.ne1 <- glmer(vote2012 ~ party2012 + age + female + 
#                            percrelchil + multi_unit_bin + med_rent_threehund +
#                            perblackbin +  perforeignbornbin + dist_pol_mi +  popdens_bin + perbachelorsbin + 
#                            (1 | geoid_new), family = binomial("logit"), 
#                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.noparty2011)

summary(fit.full2012.ne1) 

## Print to file: Table A8 (odds ratio)
## Note: Removed standard errors, stars, and "Constant" from Table A8 in the appendix
stargazer(fit.basic2012.ne1, fit.controls2012.ne1, fit.full2012.ne1,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "html", apply.coef = exp,
          out = "TableA8.html")

## Print to console: Table A8 (odds ratio)

stargazer(fit.basic2012.ne1, fit.controls2012.ne1, fit.full2012.ne1,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "text", apply.coef = exp)

#############################################
## Table A9: Odds ratios for association between no party 
## on respondent’s block in 2012 and turnout (Quasi-experimental test two)
#############################################

############################################
## Quasi-experiment 2
## Note: regression models from Figure 6
#############################################

#fit.basic2012.ne2 <- glmer(vote2012 ~ noparty2012 + 
#                             (1 | geoid_new), family = binomial("logit"), 
#                           control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.party2011)

summary(fit.basic2012.ne2) 



#fit.controls2012.ne2 <- glmer(vote2012 ~ noparty2012 + age + female +  
#                                (1 | geoid_new), family = binomial("logit"), 
#                              control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.party2011)

summary(fit.controls2012.ne2) 


#fit.full2012.ne2 <- glmer(vote2012 ~ noparty2012 + age + female +  
#                            percrelchil + multi_unit_bin + med_rent_threehund +
#                            perblackbin +  perforeignbornbin + dist_pol_mi +  popdens_bin + perbachelorsbin +
#                            (1 | geoid_new), family = binomial("logit"), 
#                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.party2011)

summary(fit.full2012.ne2) 



## Print to file: Table A9 (odds ratio)
## Note: Removed standard errors, stars, and "Constant" from Table A9 in the appendix
stargazer(fit.basic2012.ne2, fit.controls2012.ne2, fit.full2012.ne2,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "html", apply.coef = exp,
          out = "TableA9.html")

## Print to console: Table A9 (odds ratio)
stargazer(fit.basic2012.ne2, fit.controls2012.ne2, fit.full2012.ne2,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors degree (binary)",
                               "Constant"), type = "text", apply.coef = exp)

############################################
## Table A10: Association between party on 
## respondent’s block in 2012 and turnout (full sample, continuous controls)
#############################################

## version of full model with continuous controls
fit.full2012.nobin <- glmer(vote2012 ~ party2012 + age + female +  
                              percrelchil + multi_unit + med_rent_threehund + 
                              perblack +  perforeignborn + dist_pol_mi + popden + perbachelors + 
                              (1 | geoid_new), family = binomial("logit"), 
                            control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.full2012.nobin) 


## Print to file: Table A10
stargazer(fit.full2012.nobin,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (continous)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (continuous)",
                               "Percent foreign born (continous)",
                               "Distance from polling place (continuous)",
                               "Population density (continuous)",
                               "Percent bachelors degree (continuous",
                               "Constant"), type = "html", 
          out = "TableA10.html")

## Print in console: Table A10
stargazer(fit.full2012.nobin,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (continous)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (continuous)",
                               "Percent foreign born (continous)",
                               "Distance from polling place (continuous)",
                               "Population density (continuous)",
                               "Percent bachelors degree (continuous",
                               "Constant"), type = "text")

## Note: Census tract value drawn from model fit "groups": 372

############################################
## Table A11: Association between total parties on respondent’s 
## block in 2012 and turnout (full sample)
############################################

############################################
## IV -- Total parties on block in 2012
#############################################

fit.basic2012.sum <- glmer(vote2012 ~ party2012_sumall + 
                         (1 | geoid_new), family = binomial("logit"), 
                       control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.basic2012.sum) 



fit.controls2012.sum <- glmer(vote2012 ~ party2012_sumall + age + female + 
                            (1 | geoid_new), family = binomial("logit"), 
                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.controls2012.sum) 


fit.full2012.sum <- glmer(vote2012 ~ party2012_sumall + age + female +  
                        percrelchil + multi_unit_bin + med_rent_threehund + 
                        perblackbin +  perforeignbornbin + dist_pol_mi + popdens_bin + perbachelorsbin + 
                        (1 | geoid_new), family = binomial("logit"), 
                      control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.full2012.sum) 


## Print to file: Table A11
stargazer(fit.basic2012.sum, fit.controls2012.sum, fit.full2012.sum,
          covariate.labels = c("Number of parties on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent with bachelor's degree (binary)",
                               "Constant"), type = "html", 
          out = "TableA11.html")

stargazer(fit.basic2012.sum, fit.controls2012.sum, fit.full2012.sum,
          covariate.labels = c("Number of parties on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent with bachelor's degree (binary)",
                               "Constant"), type = "text")

## Note: Census tract values drawn from model fit "groups": 381, 377, 372

############################################
## Table A12: Association between total parties on respondent’s 
## (sequential dates coded as one party) block in 2012 and turnout (full sample)
############################################

############################################
## IV -- Total parties on block in 2012
## Multi-day parties coded as one 
#############################################

fit.basic2012.sum2 <- glmer(vote2012 ~ party2012_sumdate +
                         (1 | geoid_new), family = binomial("logit"), 
                       control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.basic2012.sum2) 



fit.controls2012.sum2 <- glmer(vote2012 ~ party2012_sumdate + age + female +  
                            (1 | geoid_new), family = binomial("logit"), 
                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.controls2012.sum2) 


fit.full2012.sum2 <- glmer(vote2012 ~ party2012_sumdate + age + female +  
                        percrelchil + multi_unit_bin + med_rent_threehund + 
                        perblackbin +  perforeignbornbin + dist_pol_mi + popdens_bin + perbachelorsbin + 
                        (1 | geoid_new), family = binomial("logit"), 
                      control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.full2012.sum2) 


## Print to file: Table A12
stargazer(fit.basic2012.sum2, fit.controls2012.sum2, fit.full2012.sum2,
          covariate.labels = c("Number of parties on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent with bachelor's degree (binary)",
                               "Constant"), type = "html", 
          out = "TableA12.html")

## Print to console: Table A12

stargazer(fit.basic2012.sum2, fit.controls2012.sum2, fit.full2012.sum2,
          covariate.labels = c("Number of parties on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary, cut of at median)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent with bachelor's degree (binary)",
                               "Constant"), type = "text")

## Note: Census tract values drawn from model fit "groups": 381, 377, 372

############################################
## Table A13: Association between party on respondent’s block in 
## 2012 and turnout by percent African American and foreign born in census tract
############################################

## indicator for high per black: cut-off at median value for subset > 32.43% black
## produced in code for Table 3
##data$perblackmedian <- ifelse(data$perblack >= 32, 1, 0)
##table(data$perblackmedian)

## create indicator for high per black and foreign born
data$highperblackimm <- ifelse(data$perblackmedian == 1 & data$perforeignbornbin == 1 , 1, 0)
table(data$highperblackimm)

############################################
## Regression models:
## percent black above median level
## percent foreign-born above median
#############################################


fit.full2012.highperblackimm <- glmer(vote2012 ~ party2012 + age + female +  
                                     percrelchil + multi_unit_bin + med_rent_threehund + 
                                     dist_pol_mi + popdens_bin + perbachelorsbin + 
                                     (1 | geoid_new), family = binomial("logit"), 
                                   control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data[data$highperblackimm == 1,])
summary(fit.full2012.highperblackimm) 


############################################
## Regression models:
## percent black below median level
## percent foreign-born below median
#############################################

fit.full2012.lowperblackimm <- glmer(vote2012 ~ party2012 + age + female + 
                                    percrelchil + multi_unit_bin + med_rent_threehund + 
                                    dist_pol_mi + popdens_bin + perbachelorsbin +  
                                    (1 | geoid_new), family = binomial("logit"), 
                                  control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data[data$highperblackimm ==0,])

summary(fit.full2012.lowperblackimm) 


############################################
## Interactive model
#############################################


fit.full2012.intimm <- glmer(vote2012 ~ party2012 + age + female +  
                            percrelchil + multi_unit_bin + med_rent_threehund + 
                            dist_pol_mi + popdens_bin + perbachelorsbin +
                            highperblackimm*party2012 +   
                            (1 | geoid_new), family = binomial("logit"), 
                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.full2012.intimm) 

## Print to file: Table A13
stargazer(fit.full2012.highperblackimm, fit.full2012.lowperblackimm, fit.full2012.intimm,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               #"Percent black (binary, cut of at median)",
                               #"Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors (binary)",
                               "High percent black and immigrant",
                               "At least one party on R's block in 2012 x High percent black and immigrant"
          ), type = "html", 
          out = "TableA13.html")

## Print to console: Table A13
stargazer(fit.full2012.highperblackimm, fit.full2012.lowperblackimm, fit.full2012.intimm,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               #"Percent black (binary, cut of at median)",
                               #"Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelors (binary)",
                               "High percent black and immigrant",
                               "At least one party on R's block in 2012 x High percent black and immigrant"
          ), type = "text")

## Note: Census tract values drawn from model fit "groups": 63, 309, 372

############################################
## Table A14: Association between party on respondent’s block 
## in 2012 and turnout by percentage of residents with a bachelor’s degrees in census tract
############################################

## high perc bachelors
table(data$perbachelorsbin)

############################################
## Regression model
## percent bachelor's above median level
#############################################

fit.full2012.highperbach <- glmer(vote2012 ~ party2012 + age + female + 
                                    percrelchil + multi_unit_bin + med_rent_threehund + perblackbin +
                                    perforeignbornbin + dist_pol_mi + popdens_bin + 
                                    (1 | geoid_new), family = binomial("logit"), 
                                  control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data[data$perbachelorsbin == 1,])

summary(fit.full2012.highperbach)


############################################
## Regression model
## percent bachelor's below median level
#############################################

fit.full2012.lowperbach <- glmer(vote2012 ~ party2012 + age + female +  
                                   percrelchil + multi_unit_bin + med_rent_threehund + perblackbin + 
                                   perforeignbornbin + dist_pol_mi + popdens_bin + 
                                   (1 | geoid_new), family = binomial("logit"), 
                                 control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data[data$perbachelorsbin == 0,])

summary(fit.full2012.lowperbach) 

#########################################
## Regression model
## interactive model
#############################################

fit.full2012.int <- glmer(vote2012 ~ party2012 + age + female +  
                            percrelchil + multi_unit_bin + med_rent_threehund + 
                            perblackbin +  perforeignbornbin + dist_pol_mi + popdens_bin +
                            perbachelorsbin*party2012 +   
                            (1 | geoid_new), family = binomial("logit"), 
                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.full2012.int) 

## Print to file: Table A14
stargazer(fit.full2012.highperbach, fit.full2012.lowperbach, fit.full2012.int,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelor's degree (binary)",
                               "At lest one party on R's block in 2012 x Percent bachelor's",
                               "Constant"
                               
          ), type = "html", 
          out = "TableA14.html")

## Print to console: Table A14
stargazer(fit.full2012.highperbach, fit.full2012.lowperbach, fit.full2012.int,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Median rent (continuous; unit = $300)",
                               "Percent black (binary)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelor's degree (binary)",
                               "At lest one party on R's block in 2012 x Percent bachelor's",
                               "Constant"
                               
          ), type = "text")

## Note: Census tract values drawn from model fit "groups": 193, 179, 372


############################################
## Table A15: Association between party on respondent’s block 
## in 2012 and turnout by median rent in census tract
############################################

## cut-off at median value for subset > 858
summary(data$rent_bin) # 80 missing
## remove nas
data.3 <- data[!is.na(data$rent_bin),] ## 99291 obs
table(data.3$rent_bin)

############################################
## Regression model
## Avg. rent above median level
#############################################


fit.full2012.highrent <- glmer(vote2012 ~ party2012 + age + female +  
                                 percrelchil + multi_unit_bin + perblackbin + 
                                 perforeignbornbin + dist_pol_mi + popdens_bin + perbachelorsbin +  
                                 (1 | geoid_new), family = binomial("logit"), 
                               control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.3[data.3$rent_bin == 1,])

summary(fit.full2012.highrent) 

############################################
## Regression model
## Avg. rent below median level
#############################################


fit.full2012.lowrent <- glmer(vote2012 ~ party2012 + age + female +  
                                percrelchil + multi_unit_bin + perblackbin +
                                perforeignbornbin + dist_pol_mi + popdens_bin + perbachelorsbin +
                                (1 | geoid_new), family = binomial("logit"), 
                              control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data.3[data.3$rent_bin == 0,])

summary(fit.full2012.lowrent) 

############################################
## Regression model
## interactive model
#############################################

fit.full2012.int <- glmer(vote2012 ~ party2012 + age + female +  
                            percrelchil + multi_unit_bin + 
                            perblackbin +  perforeignbornbin + dist_pol_mi + popdens_bin + perbachelorsbin +
                            rent_bin*party2012 +   
                            (1 | geoid_new), family = binomial("logit"), 
                          control = glmerControl(optimizer = "nloptwrap", calc.derivs = FALSE), data = data)

summary(fit.full2012.int) 

## Print to file: Table A15
stargazer(fit.full2012.highrent, fit.full2012.lowrent, fit.full2012.int,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Percent black (binary)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelor's degree (binary)",
                               "Median rent (binary)",
                               "At lest one party on R's block in 2012 x High median rent",
                               "Constant"
                               
          ), type = "html", 
          out = "TableA15.html")

## Print to console: Table A15
stargazer(fit.full2012.highrent, fit.full2012.lowrent, fit.full2012.int,
          covariate.labels = c("At least one party on R’s block in 2012",
                               "Age",
                               "Female",
                               "Percent related child (continuous)",
                               "Percent multi-unit (binary)",
                               "Percent black (binary)",
                               "Percent foreign born (binary)",
                               "Distance from polling place (continuous)",
                               "Population density (binary)",
                               "Percent bachelor's degree (binary)",
                               "Median rent (binary)",
                               "At lest one party on R's block in 2012 x High median rent",
                               "Constant"
                               
          ), type = "text")
## Note: Census tract values drawn from model fit "groups": 189, 183, 372
